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Abstract 

We analyze low field hysteresis close to the demagnetized state in disordered 
ferromagnets using the zero temperature random-field Ising model. We solve 
the demagnetization process exactly in one dimension and derive the Rayleigh 
law of hysteresis. The initial susceptibility a and the hysteretic coefficient b 
display a peak as a function of the disorder width. This behavior is confirmed 
by numerical simulations d = 2, 3 showing that in limit of weak disorder 
demagnetization is not possible and the Rayleigh law is not defined. These 
results are in agreement with experimental observations on nanocrystalline 
magnetic materials. 

PACS numbers: 75.60.Ej, 75.60.Ch, 64.60.Ht, 68.35.Ct 
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I. INTRODUCTION 



Ferromagnetic materials display hysteresis under the action of an external field and the 
magnetization depends in a complex way on the field history. In order to define magnetic 
properties unambiguously, it is customary to first demagnetize the material, bringing it to a 
state of zero magnetization at zero field. This can be done, in practice, by the application 
of a slowly varying AC field with decreasing amplitude. In this way, the system explores a 
complex energy landscape, due to the interplay between structural disorder and interactions, 
until it is trapped into a low energy minimum. This demagnetized state is then used as a 
reference frame to characterize the magnetic properties of the material. 

The hysteresis properties at low fields, starting from the demagnetized state, have been 
investigated already in 1887 by Lord Rayleigh [p]], who found that the branches of the 
hysteresis loop are well described by parabolas. In particular, when the field is cycled 
between ±H* , the magnetization M follows M = (a + bH*)H ± b{{H*) 2 - H 2 )/2, where 
the signs ± distinguish the upper and lower branch of the loop. Consequently the area of 
the loop scales with the peak field H* as W = A/3b(H*) 3 and the response to a small field 
change, starting from the demagnetized state is given by M* = a(H*) ± b(H*) 2 0. 

The Rayleigh law has been widely observed in ferromagnetic materials ||, but also in 
ferroelectric ceramics ^M- The current theoretical interpretation of this law is based on 
a 1942 paper by Neel H, who derived the law formulating the magnetization process as 
the dynamics of a point (i.e. the position of a domain wall) in a random potential. In 
this framework, the initial susceptibility a is associated to reversible motions inside one 
of the many minima of the random potential, while the hysteretic coefficient b is due to 
irreversible jumps between different valleys. Successive developments and improvements 
have been devoted to establish precise links between Neel random potential and the material 
microstructure [§-0, but in several cases the issue is still unsettled. For instance, the initial 
permeability of nano crystalline materials typically displays a peak as a function of the grain 
size |T(J, heat treatment JXTJJT^] or alloy composition J10|JT^]. This behavior can be associated 
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to changes in the disordered microstructure, but can not be accounted for by Neel theory 
that predicts a monotonic dependence of a on the disorder width |J. 

The zero temperature random-field Ising model (RFIM) has been recently used to de- 
scribe the competition between quenched disorder and exchange interactions and their ef- 
fect on the hysteresis loop |14j| . In three and higher dimensions, the model shows a phase 
transition between a continuous cycle for strong disorder and a discontinuous loop, with a 
macroscopic jump, at low disorder. The two phases are separated by a second order crit- 



ical point, characterized by universal scaling laws |]14HT6[| . A behavior of this kind is not 
restricted to the RFIM but has also been observed in other models, with random bonds or 
random anisotropics [|17| and vectorial spins ||18|| . In addition, a similar disorder induced 



phase transition in the hysteresis loop has been experimentally reported for a Co-CoO bi- 
layer [fL9). Thus the RFIM provides a tractable model for a more generic behavior: the 



model has been solved exactly in one dimension pO| , |2T| and on the Bethe lattice p2| , p3| , 
while mean-field theory [14| and renormalization group [T^| have been used to analyze the 
transition. 

Here, we use the RFIM to analyze the demagnetization process and investigate the 
properties of the hysteresis loop at low fields. Along the lines of Refs. [^T],^3J, we compute 



the demagnetization cycles exactly in one dimension and derive the Rayleigh law, obtaining 
a and b as function of disorder and exchange energies. Next, we analyze the problem 
numerically in higher dimensions (i.e. d = 2 and d = 3) where exact results are at present 



not available. In d = 3, we find that the disorder induced transition JL4J], defined on the 
saturation loop, is also reflected by the Rayleigh loops: in the weak disorder phase the 
system can not be demagnetized, as the final magnetization coincides with the saturation 



magnetization. A similar behavior has been recently obtained analyzing subloops P^| . In 
the high disorder phase, however, a demagnetization process is possible and hysteresis loops 
are still described by the Rayleigh law. Above the transition, the dependence of a and b on 
disorder is qualitatively similar in all dimensions, displaying a peak and decreasing to zero 
for very strong disorder in agreement with experiments [|10H T3||. 



II. THE RANDOM FIELD ISING MODEL 



In the RFIM, a spin Sj = ±1 is assigned to each site i of a d— dimensional lattice. The 
spins are coupled to their nearest-neighbors spins by a ferromagnetic interaction of strength 
J and to the external field H. In addition, to each site of the lattice it is associated a random 
field hi taken from given probability distribution p(h). In the following we will mainly focus 
on a Gaussian with variance R (i.e. p(h) = exp(—h 2 /2R 2 )/^/2nR), but we will also consider 
a rectangular distribution. The Hamiltonian thus reads 

H = - y £Js i s j - y E(H + h i )s i , (1) 

where the first sum is restricted to nearest-neighbors pairs. The dynamics proposed in 



Ref. [25] and used in Refs. |14|-|T6] is such that the spins align with the local field 



Si =sign{JY^Sj + hi + H). (2) 

3 

In d = 1, a spin with n neighbors up (n = 0,1,2), will be up at the field H with 
probability: 

r+oo 

Pn{H) = / p(hi)dhi (3) 

J2(l-n)J-H 

When a spin flips up the local field of its neighbors is raised by 2J so that it can happen 
that one or both of the two neighbors flip up. In this way a single spin flip can lead the 
neighboring spins to flip, eventually triggering an avalanche. 



It has been shown that the RFIM obeys return-point memory |L4| : if the field is increased 
adiabatically the magnetization only depends on the state in which the field was last reversed. 
This property has been exploited in d — 1 and in the Bethe lattice to obtain exactly the 



saturation cycle and the first minor loops |21j|. In the next section we will briefly recall the 



results reported in Ref. |^TJ and we will then proceed with a general derivation for nested 
minor loops. 
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III. SATURATION LOOP AND FIRST RETURN CURVES 



To obtain the saturation loop, we start from the initial condition Sj = — 1 at H = — oo 
and we will raise the field up to Hq. We are thus moving on the lower half of the major 



hysteresis loop. Following Ref. |2T|, we define the conditional probability Uq that a spin 
flips up at Hq before a given nearest neighbor. To compute Uq, we take advantage of the 
translational invariance of the system. There are only two ways to flip up a spin in i keeping 
the spin in i — 1 down. The two contributions yield U = pi(H )U + p (H )[l — U ], from 
which we obtain 

= Po{H ) 
1 - [ Pl (H ) - p (H )} {) 

The probability that a spin is up at field Ho is 

p(H ) = U 2 p 2 (H ) + 2U (1 - U ) Pl (H ) + (1 - U ) 2 Po (H ) (5) 

and the magnetization per spin M(H ) is simply M(H Q ) = 2p(H ) — 1. In Fig.|l] we show 
the saturation loop for a Gaussian distribution of random fields. 

If the field is reversed from a finite value H , we have a new situation and the system 
departs from the saturation curve. It is possible to show that if the field changes from Hq 
to H\ = Hq — 2 J the magnetization reaches the upper saturation loop again. Thus we can 
restrict the analysis to fields included in [Hq — 2 J, Hq}. The first return curve can be obtained 
counting the spins that were up at H and are down at Hi. To this end, we introduce Di as 
the conditional probability that a spin is down if its neighbor is up. Following similar steps 



as for Uq ||21|| , we obtain: 



1 - \pi(H ) - pi(Hi)] ' {) 
where f(H ) = f/ [l — pi(H )} + (1 — Uo)[l — p2(H )]. At this point it is straightforward to 
write the probability p(Hi) that a spin is up at Hi: 

p(Hi) =p(H )-U 2 ip 2 (H )-p 2 (Hi)} (7) 
+2U D 1 \ P i(Ho)-pi(H 1 )}+D 2 1 \p {Ho)- Po (H 1 )} 



which is simply related to the magnetization. 



IV. DEMAGNETIZATION 



Here, we extend the approach of Ref. []2T[ to more general field histories, treating explic- 
itly the demagnetization process: the external field is changed through a nested succession 

H = H -> #1 -> H 2 -> fln- -> 0, with F 2n > # 2 n+2 > 0, H 2n -i < H 2n +i < and 

dH = H 2n — H 2n+2 — ► 0. The initial value i?o should correspond to complete saturation, but 
we discussed above that as long as H n > J the magnetization M n = M(H n ) simply follows 
the saturation curve, so that we can set H Q = J. 

As in the previous section, the key quantity to compute is the conditional probability 
U 2n that a spin flips up before its nearest neighbor when the field is increased from H 2n _\ 
to H 2n . Similarly on the descending part of the loops we define -D 2n +i as the conditional 
probability that a spin flips down before its nearest neighbor when the field is decreased 
from H 2n to H 2n+ i. Enumerating all possible spin histories, we find recursion relations for 
the conditional probabilities which read as |p6| 

U2n-2[Pl{H 2n ) ~ p 1 (H 2n - 2 )] + £>2n-l[Po(#2n) ~ Po{H 2n - 2 )) 



U- 



2n 



2n-2 



+ 



D 2n +\ — D 2n -i + 



l-\p 1 {H 2n )-p 1 {H 2n ^ 1 )] 

£>2n-l[Pl(#2n-l) -p\{H 2n+1 )} + U 2n \p 2 {H 2n ^) -p 2 (H 2n+ i)] 



(8) 



1 - \pi(H 2n ) -pi(if 2 n+l)] 
The derivation of Eqs. @ is a little involved and we thus report it in the Appendix. 

The magnetization as a function of the peak field is given by 



M 2n = M 2n _x + 2Ul n [p 2 (H 2n 
+4[/ 2nJ D 2n _ 1 [p 1 (# 2n ; 
+2 J D 2 V 1 [p (tf 2n ; 



and a similar expression holds for M 



2n+l- 



In the limit H- 



2n-2 



H- 



2n 



dH 



0, H 



2n 



-p 2 (H 2n „i)) 

-Pl(i?2„-l)] 

-Po{H 2n -i)} 
H* and H 2n _ x 



(9) 



-H* , the recursion 



relations in Eqs. |g become a pair of differential equations |27 



r au_ 

~dH* 
db 



1 



I dH* 



1-Q 
1 



p(H*)D + p(2J - H*)U 



(10) 



p{H*)U - p{2J - H*)D , 

where Vt = f H H «p(h')dh' and D(H) = D(—H). The boundary conditions are given by the 
conditional probabilities on the saturation loop (i.e. U(J) = D(J) = U (J) = 1/2) and the 
solution reads 



p(h') + p(2J - h') 



dti 



(11) 



ih* i - n(h') 

Once the conditional probability U is known, it is straightforward to compute the magneti- 
zation as a function of the peak field H* from Eq. |9|, noting that M(—H*) = —M(H*). Inner 
loops starting from the demagnetization curve (i.e. Eq. |9|) can also be computed exactly. In 
Fig. m we report the demagnetization curve and a few inner loops for a system with Gaus- 
sian random field distribution with unit variance. The analytical results are compared with 
numerical simulations, performed on a lattice with L = 5 10 5 spins, using a single realization 
of the disorder. The perfect agreement between the curves confirms that the magnetization 
is self- averaging, as assumed throughout the calculations. 



V. RAYLEIGH LAW 

To analyze low field hysteresis we first substitute in Eq. [| H2 n and B.2 n -\ with H* and 
—H*. If we start to reverse the field from Hq = J and we cycle the field symmetrically around 
H* = 0, the process displays the symmetry M(H*) = -M(-H*) and U(H*) = D{H*). 
Thus we can reduce Eq. |] to 

M(H*) = 2U 2 (H*) £ \p k (H*) - Pk (-H% (12) 

fc=0 

Now we can expand M(H*) around H* = 0. In this limit we have 

, 2H*p(2J) if k = 0,2 
\p k (H*)- Pk (-H*)]^{ (13) 

2H*p(0) if k = 1 



and 

U 2 (H*) ~ U 2 (0)[1 + 2H*(p(0) + p{2J))\. (14) 

Collecting Eq. [L3]and Eq. [b|m Eq. [12], we obtain M ~ aH* + b(H*) 2 recovering the Rayleigh 
expression with 



a = 4f/ 2 (0)[p(0) +p(2J) 
b = 4f/ 2 (0)[p(0) + p(2J) 



(15) 



An expansion can also be performed for minor loops on the demagnetization curve (i.e 
cycling H between ±H*), yielding M = (a + bH*)H ± b((H*) 2 — H 2 )/2, which coincides 
with the Rayleigh law. 

In Fig. Qa we report the values of a and b for a Gaussian distribution of random fields as 
a function of the disorder R, showing that both components of the susceptibility display a 

maximum in R. To identify the low and strong disorder behavior of the susceptibilities, we 

2 2 

perform an asymptotic expansion and we obtain for R — > oo that a ~ and b 



2itR " ' irR 2 ' 

For R — > 0, we obtain: a ~ ( )e~^R I and 6 ~ [ == ) — e - ^. Finally in Fig. |2|b we 

1 evrjy \e7ijV2^J R U 



report a and 6 obtained with a rectangular distribution of random fields. The derivation of 
these results is reported in appendix B. 



VI. SIMULATIONS IN D = 2, 3 

Next, we turn our attention to high dimensional system, for which analytical results 
are not available. In order to obtain unambiguously the demagnetized state for a given 
realization of the disorder, one should perform a perfect demagnetization. This is done in 
practice changing the field by precisely the amount necessary to flip the first unstable spin. 
In this way, the field is cycled between —H* and H* and H* is then decreased at the next 
cycle by precisely the amount necessary to have one avalanche less than in the previous cycle. 
This corresponds to decrease H* at each cycle by an amount dH, with dH — > + . The perfect 
demagnetization algorithm allows to obtain a precise characterization of the demagnetized 
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state but it is computationally very demanding. Thus we resort to a different algorithm 
which performs an approximate demagnetization: instead of cycling the field between — H* 
and H* we just flip the field between these two values and then decrease H* by a fixed 
amount dH . We have checked that with a reasonably small dH (i.e. dH < 10~ 3 ) the 
demagnetization curve is quite insensitive to the algorithm used. 

As we discussed above, it is well established that in d = 3 the saturation loops reveal 
a phase transition at R c ~ 2.16 for J = 1 ]16| (the transition is not present in d = 1, 



while in d = 2 the issue is controversial [lq]). We find that the transition is reflected also 
in the Rayleigh loops: in Fig. |3] we report the final magnetization computed using 
the demagnetization algorithm for different values of R. For strong disorder R > R c , we 
see that ~ as expected, but as R < R c the demagnetization curve tends to the 
saturation magnetization and — > ±1. The transition becomes sharper as the system 
size is increased, indicating that demagnetization is possible only for R > R c (see also 
Ref. [p^]). We notice here that two scenarios are possible for L —>■ oo as R — > R~ . The 
first possibility is that scales continuously to zero as (R c — R)@ and the second is that 
the transition is discontinuous (i.e. — > M* > 0). The present numerical results do not 
allow to distinguish between these two cases, but a recent analysis of the RFIM on the Bethe 
lattice is in favor of the first alternative p8| . 

From the demagnetization curve, the Rayleigh parameters can be estimated plotting 
M — Moo/H vs H and fitting the linear part of the curve close to H = (see Fig. As we 
show in Fig. |] the demagnetization curve is basically independent from the system size, once 
the magnetization has been shifted by M^. Thus we expect that the Rayleigh parameters 
be also independent on L. In Fig. |^ we report the values of a and b obtained numerically in 
d = 2 and d = 3 for different values of R, using systems of sizes (L = 100) 2 and (L = 50) 3 . 
The results are qualitatively similar to those obtained exactly in d = 1: the curve displays 
a peak for intermediate disorder and decrease to zero for weak and strong disorder. 
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VII. DISCUSSION 



In this paper we have discussed the demagnetization properties of the RFIM in d = 1, 2, 3. 
In d = 1 it is possible to compute exactly the demagnetization curve and obtain an expression 
for the Rayleigh parameters. We find that a and b display a peak in the disorder R. This 
result is confirmed by numerical simulations in d = 2, 3, where analytical results are not 
available. In addition, in d — 3 the disorder induced phase transition strongly affects the 
demagnetization process: for R < R c it is not possible to demagnetize the system anymore. 

It is interesting to compare our theoretical results with experiments on nanocrystalline 
materials. It has been reported that the initial susceptibility in several cases displays a 
peak as the heat treatment or the alloy composition are varied HlQHilfl . The peak is usually 
associated to changes in the microstructure, which induce a competition between the disorder 
present in grain anisotropics and inter-grain interactions mediated by the amorphous matrix 
|T0| . Notice that a similar behavior can not be reproduced by Neel theory, where the initial 



susceptibility is decreasing with the width of the disorder potential ||. On the other hand, we 
see here that the behavior is well captured by the RFIM, that allows to analyze the the effect 
of the disorder-exchange ratio R/J. For weak disorder, we have a few large domains and 
the susceptibility is dominated by domain wall dynamics. When the disorder is increased, 
the number of domains (and domain walls) also increases and so does the susceptibility. 
Increasing the disorder further leads to a complete breakup of the domains and the response 
is dominated by single spin flips in low random-field regions with a progressive decrease of 
the susceptibility. 

A detailed understanding of the demagnetization process and low field hysteresis has 
important implications also from a purely theoretical point of view. When a disordered 
system is demagnetized, it explores a complex energy landscape until it finds a metastable 
minimum. It would be interesting to compare the statistical properties of the demagnetized 
state, with those of the ground state of the system [p5|l . The analysis of the ground state of 
disordered systems has received a wide attention in the past few years, due to the connections 
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with general optimization problems, and the RFIM is one of the typical model used to test 
ground state algorithms ||29|| . Demagnetization could provide a relatively simple way to 
obtain a low energy state that can be useful for optimization procedures. We are currently 
pursuing investigations along these lines [P0f| . 
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APPENDIX A: DERIVATION OF THE RECURSION RELATIONS 

Here, we derive recursion relations for the conditional probabilities U 2n and D 2n+ i as a 
function of the previous magnetization history. Let us first consider the case of D 2n +i- the 
field from H 2n -i reaches H 2n and is then decreased again up to H 2n+ i. The weight of the 
fraction of spins that at field H 2n +i flip down before their neighbor is given by 

D 2n +i = D 2n -l ~ C,2n + C2n+1, (Al) 

where ( 2n is the weight of the fraction of spins that were down at H 2n _i before a fixed nearest 
neighbor and flip up at H 2n , while ( 2n +i is the weight of the fraction of spins contributing 
to C, 2n which flip again down at H 2n +i- 

To compute ( 2n , we consider the spins that at the field H 2n -i are down before their 
neighbor (for instance, we can say that the spin z-th is down before the spin in site % — 1) 
and are up at the field H 2n . Since we fixed up the spin in site i — 1, the spin in site i + 1 
can be either up or down. If the spin in i + 1 is up when the spin i flips up, it contribution 
to ( 2n with 

U 2n [p 2 (H 2n ) -p 2 (#2n-l)]- 



11 



If the spin in site i + 1 is down when the spin i flips up, we obtain 

D 2n -i[pi(h 2n ) -pi(h 2n -i)}- 

Indeed, \p n (H 2n ) — p n (H 2n -i)} is the probability that a spin with n up nearest neighbors is 
up at H 2n but not at i^n-i, while D 2n _i,U 2n are respectively the conditional probabilities 
that the spin in site i + 1 is down or up if the spin in site i is down. Adding the two 
contributions, we obtain 

= {D 2n -i\Pi{H 2n ) -pxiH^)) + U 2n \p 2 (H 2n ) -p 2 (#2n-i)]}- (A2) 

The derivation of C2n+i follows similar steps: we count the spins that are up at H 2n and are 
again down at H 2n+1 . If the spin in the site i + 1 is up at H 2n+ i, the spin in i is up at H 2n 
and is down at H 2n+ i with probability 

U 2n [p 2 (H 2n ) -p 2 (H 2n+1 )]. 

Finally, we analyze the case in which the spin in site i + 1 is already down when the spin i 
flips down. The weight of this configuration is 

D 2n+ i[pi{h 2n ) -Pi{h 2n+1 )], 

so that ( 2n+ i is given by 

(2n+i = {u 2n [p2(H 2n ) - p 2 (H 2n+1 )} + D 2n+1 [ Pl (H 2n ) - Pl (H 2n+1 )}}. (A3) 

Substitute these two expressions in Eq. ( |Al| ) we obtain the second of Eqs. ||. We can then 
derive a similar equation for U 2n (First of Eqs. |8|) following the same method as the one 
employed above to calculate D 2n+ i. 

APPENDIX B: THE CASE OF THE RECTANGULAR DISTRIBUTION 

It is also instructive to consider the case of a rectangular distribution of random fields 
(i.e. p(x) = 1/2 A if A < x and zero otherwise), since all the calculations can be carried 
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out explicitely. As usual, we cycle the field around H = and we take Hq = J. The 
calculation should be divided in several cases, depending on the value of A. 

(i) For A > 3 J, we have p(x) = p(2J - x) = 1/2A, so that p k (H*) -p k (-H*) = H* / A 
and Eq. ([□]), reduces to 



U 2 (H* 



A- J 
4 \A-H* 



:bi) 



Inserting these results in Eq. (|12"1), we obtain: 

M(H*" ' A 1 



A-H* 



H 
A' 



(B2) 



Expanding Eq. B3, we obtain the values for a and b 



a =A 



1 - ± 

1 A 



(B3) 



b — 1 — 



1 



(ii) For 2 J < A < 3 J, U 2 (0) is still given by Eq. QBT]) but differs from the 

previous case. The magnetization is now given by 

2 3H *~ A l J+A H*>A-2J 

lA (B4) 



M(H*) = (^y) f H* < A - 2J. 

The expansion around if* is thus still given by Eq. (B3). 

(iii) The behavior for J < A < 2 J is again different: close to H* = the peak magneti 
zation is not given by Eq. ([12]), but for H* < 2 J — A can be written as 

(A - JfH* 



M(H* 



so that expanding we obtain 



AJA(A-H*)' 



(A- J) 2 



(B5) 



4A 2 J 

(A-J) 2 
2A3J 



(B6) 



(iv) Finally for A < J there is no hysteresis and thus the Rayleigh law is not defined. 
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FIG. 1. Exact expressions for the saturation cycle (thin lines), the demagnetization curve (thick 
lines) and a few minor loops (dotted lines) for J = 1 and R = 1. The points are the results of a 
numerical simulation with L = 5 10 5 spins and a single realization of the disorder. 
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FIG. 2. The reversible susceptibility a and the hysteretic coefficient b computed exactly in 
d = 1 for (a) a Gaussian distribution of random fields and (b) a rectangular distribution. 
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FIG. 3. The absolute value of the final magnetization \M 00 \ as a function of R, obtained from 
numerical simulations in d = 3. For strong disorder |Moo| = as expected, while for weak disorder 
the final magnetization coincides with the saturation value. The transition between the two types 
of behavior becomes sharper as the system size is increased. 
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FIG. 4. The demagnetization curve can be used to obtain an estimate of the Rayleigh param- 
eters. Notice the absence of system size dependence. These results are obtained in d = 3. 
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FIG. 5. The reversible susceptibility a computed exactly in d = 1 is compared with numerical 
results in d = 2 and d = 3. In the inset we show a similar plot for the parameter b. 
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